home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
SGI Developer Toolbox 6.1
/
SGI Developer Toolbox 6.1 - Disc 4.iso
/
public
/
bit
/
src
/
iprocess.c
< prev
next >
Wrap
C/C++ Source or Header
|
1994-08-01
|
23KB
|
905 lines
/*
* $Id: iprocess.c,v 0.91 1994/02/20 00:53:18 zhao Pre-Release $
*
*. This file is part of BIT shareware package. After the two weeks of
* free evaluation period, you are encouraged (required) to register
* your copy for a small registration fee, which is $35 for personal use
* and $50 for commercial, government and institutional use.
*
* Copyright(c) 1993, 1994 by T.C. Zhao.
* All rights reserved.
*
* Permission to use, copy, and distribute this software in its entirety
* for non-commercial purposes is hereby granted, provided that the
* above shareware and copyright notices and this permission notice
* appear in all copies and their documentation.
*
* This software may be modified for your own use, but modified versions
* may not be distributed without prior consent of the author.
*
* This software is provided "as is" without expressed or implied
* warranty of any kind.
*
*.
*
* Simple image processing routines. Not very isolated, but there
* are no-interactive tasks to be performed by any of these routines.
*/
#if !defined(lint) && defined(F_ID)
char *id_ip = "$Id: iprocess.c,v 0.91 1994/02/20 00:53:18 zhao Pre-Release $";
#endif
#include "bit.h"
/********************************************************************
* General convlution routines.
* could be further optimized/hacked to go faster, but undoubtly
* will be uglier.
*******************************************************************/
#define vectorP3(cm, m, c) \
(cm[0]*m[c-1] + cm[1] * m[c] + cm[2]*m [c + 1]) \
#define vectorP5(cm, m, c) \
(cm[0]*m[c-2]+cm[1]*m[c-1]+cm[2]*m[c]+ \
cm[3]*m[c+1]+cm[4]*m[c+2])
/* Perform a 3x3 convolution at m(r,c) */
#define conv3x3(cm,m,r,c) \
(vectorP3(cm[0], m[r-1], c) + \
vectorP3(cm[1], m[r], c) + \
vectorP3(cm[2], m[r + 1], c))
/* may require compiler flag -Wf,-XNh2000 to */
#if 0
#define conv5x5(cm,m,r,c) \
(vectorP5(cm[0], m[r - 2], c) + vectorP5(cm[1], m[r - 1], c) + \
vectorP5(cm[2], m[r] , c) + vectorP5(cm[3], m[r + 1], c) + \
vectorP5(cm[4], m[r + 2], c))
#else
static int
conv5x5(register int **cm, register pc_t **m,
register int r, register int c)
{
return (vectorP5(cm[0], m[r - 2], c) +
vectorP5(cm[1], m[r - 1], c) +
vectorP5(cm[2], m[r], c) +
vectorP5(cm[3], m[r + 1], c) +
vectorP5(cm[4], m[r + 2], c));
}
#endif
#define PC_check(pc) \
do { \
if(pc < 0) pc = 0; \
else if((pc /= weight) > PCMAXV) \
pc = PCMAXV; \
} while (ZERO)
/*** thresholding ****/
#define TRSHD(new, old, t, i) \
((t && t[i]) ? ((Abs((int)(new)-(int)(old)) >= t[i]) ? \
(new):(old)) : (new))
/****************************************************************
* Given RGB matrices and the convolution kernel, do the
* convolution.
****************************************************************/
static int
rgb_convolve(register pc_t **rm, register pc_t **gm, register pc_t **bm,
int nrow, int ncol, register int **mat, int cm_col, int cm_row,
register int *threshold, const char *what)
{
register long weight = 0;
register int rnew, gnew, bnew;
register int row, col;
int i, j;
long rlines;
/* get normalization factor */
for (row = 0; row < cm_row; row++)
for (col = 0; col < cm_col; col++)
weight += mat[row][col];
if (weight <= 0)
{
Bark("Convolve", "Bad matrix weight: %ld", weight);
return -1;
}
rlines = progress_report(what, nrow - cm_row / 2);
if (cm_col == 3 && cm_row == 3)
{
nrow--;
ncol--;
for (row = 1; row < nrow; row++)
{
REPORT(row, rlines);
for (col = 1; col < ncol; col++)
{
rnew = conv3x3(mat, rm, row, col);
gnew = conv3x3(mat, gm, row, col);
bnew = conv3x3(mat, bm, row, col);
PC_check(rnew);
PC_check(gnew);
PC_check(bnew);
rm[row][col] = TRSHD(rnew, rm[row][col], threshold, 0);
gm[row][col] = TRSHD(gnew, gm[row][col], threshold, 1);
bm[row][col] = TRSHD(bnew, bm[row][col], threshold, 2);
}
}
}
else if (cm_col == 5 && cm_row == 5)
{
nrow -= 2;
ncol -= 2;
for (row = 2; row < nrow; row++)
{
REPORT(row, rlines);
for (col = 2; col < ncol; col++)
{
rnew = conv5x5(mat, rm, row, col);
gnew = conv5x5(mat, gm, row, col);
bnew = conv5x5(mat, bm, row, col);
PC_check(rnew);
PC_check(gnew);
PC_check(bnew);
rm[row][col] = TRSHD(rnew, rm[row][col], threshold, 0);
gm[row][col] = TRSHD(gnew, gm[row][col], threshold, 1);
bm[row][col] = TRSHD(bnew, bm[row][col], threshold, 2);
}
}
}
else
{
int ccol, ii, jj;
int cm_2_row = cm_row / 2;
int cm_2_col = cm_col / 2;
for (row = cm_2_row; row < nrow - cm_2_row; row++)
{
REPORT(row, rlines);
for (col = cm_2_col; col < ncol - cm_2_col; col++)
{
rnew = gnew = bnew = 0;
ccol = col - cm_2_col;
for (i = 0; i < cm_row; i++)
{
ii = row - cm_2_row + i;
for (jj = ccol, j = 0; j < cm_col; jj++, j++)
{
rnew += mat[i][j] * rm[ii][jj];
gnew += mat[i][j] * gm[ii][jj];
bnew += mat[i][j] * bm[ii][jj];
}
}
/* clamp to minmax */
PC_check(rnew);
PC_check(gnew);
PC_check(bnew);
rm[row][col] = TRSHD(rnew, rm[row][col], threshold, 0);
gm[row][col] = TRSHD(gnew, gm[row][col], threshold, 1);
bm[row][col] = TRSHD(bnew, bm[row][col], threshold, 2);
}
}
}
remove_progress_report();
return 0;
}
/*****************************************************************
* Similar to rgb_convolve except that the input is in grayscale.
* Note gm & bm are un-used
******************************************************************/
/* ARGSUSED */
static int
gray_convolve(register pc_t **rm, pc_t **gm, pc_t **bm,
int nrow, int ncol, register int **mat, int cm_col, int cm_row,
register int *threshold, const char *what)
{
register long weight = 0;
register int rnew;
register int row, col, i, j, ii, jj;
long rlines;
/* get normalization factor */
for (row = 0; row < cm_row; row++)
for (col = 0; col < cm_col; col++)
weight += mat[row][col];
if (weight <= 0)
{
Bark("Convolve", "Bad matrix weight: %ld", weight);
return -1;
}
rlines = progress_report(what, nrow - cm_row / 2);
if (cm_col == 3 && cm_row == 3)
{
nrow--;
ncol--;
for (row = 1; row < nrow; row++)
{
REPORT(row, rlines);
for (col = 1; col < ncol; col++)
{
rnew = conv3x3(mat, rm, row, col);
PC_check(rnew);
rm[row][col] = TRSHD(rnew, rm[row][col], threshold, 0);
}
}
}
else if (cm_col == 5 && cm_row == 5)
{
nrow -= 2;
ncol -= 2;
for (row = 2; row < nrow; row++)
{
REPORT(row, rlines);
for (col = 2; col < ncol; col++)
{
rnew = conv5x5(mat, rm, row, col);
PC_check(rnew);
rm[row][col] = TRSHD(rnew, rm[row][col], threshold, 0);
}
}
}
else
{
int ccol;
int cm_2_row = cm_row / 2;
int cm_2_col = cm_col / 2;
for (row = cm_2_row; row < nrow - cm_2_row; row++)
{
REPORT(row, rlines);
for (col = cm_2_col; col < ncol - cm_2_col; col++)
{
rnew = 0;
ccol = col - cm_2_col;
for (i = 0; i < cm_row; i++)
{
ii = row - cm_2_row + i;
for (jj = ccol, j = 0; j < cm_col; jj++, j++)
{
rnew += mat[i][j] * rm[ii][jj];
}
}
/* clamp to minmax */
PC_check(rnew);
rm[row][col] = TRSHD(rnew, rm[row][col], threshold, 0);
}
}
}
remove_progress_report();
return 0;
}
/*************************************************************
* Global interface to convolution routines
*
* do the convolution on an image or a part of it.
*************************************************************/
int
img_convolv(IPTR im, register int **m, int nrow, int ncol,
int *t, const Rect_t * sr, const char *msg)
{
register pc_t **rm, **gm, **bm;
int status;
int total = sr->h * sr->w;
rgba_t **subi;
/* always demand proper region */
if (sr->w < ncol || sr->h < nrow)
return -1;
fl_deactivate_all_forms();
/*
* smoothing can only be done to RGB images. Turn maped gray image to
* gray and all others to RGB
*/
if (!IS_CPACK(im))
{
(void) img_convert_type(im, IS_GRAY(im) ? T_GRAY : T_RGBA);
/* re-display */
im->io->display(im, -1, 0);
}
/* check if to be done on the entire image */
if (!equal_rect(sr, img_rect(im)))
{
/* not whole image */
subi = get_subimage(im, sr->x, sr->y, sr->w, sr->h);
/* gm & bm are un-used n convolution is grayscale image */
rm = get_mat(sr->h, sr->w, sizeof(pc_t));
gm = get_mat(sr->h, sr->w, sizeof(pc_t));
bm = get_mat(sr->h, sr->w, sizeof(pc_t));
/* unpack RGBA to R, G, B. */
unpack_mat(subi[0], rm[0], gm[0], bm[0], total);
}
else
{
subi = im->mraster;
img_get_rgb(im);
rm = get_RM(im);
gm = get_GM(im);
bm = get_BM(im);
}
status = (IS_GRAY(im) ? gray_convolve : rgb_convolve)
(rm, gm, bm, sr->h, sr->w, m, nrow, ncol, t, msg);
/* pack_mat detects gray & rgba by checking if R == G */
pack_mat(subi[0], rm[0], IS_GRAY(im) ? rm[0] : gm[0], bm[0], total);
/* change the image in core */
put_subimage(im, subi, sr, 1);
/* change the image in framebuffer */
set_current_window(win_id);
Rectwrite(sr->x, sr->y, sr->x + sr->w - 1, sr->y + sr->h - 1, subi[0]);
/* do the clean up */
if (subi == im->mraster)
{
img_free_rgbmem(im);
}
else
{
free_mat(rm);
free_mat(gm);
free_mat(bm);
free_mat(subi);
}
fl_activate_all_forms();
return status;
}
/*
* sum of two derivatives and convert to colormapped image only looks ok
* followed by a thresholding with coledit
*/
#define get_derivative(r0, r1, r2, i) \
do { \
d1 = (int) r0[i + 1] - (int) r0[i - 1] + \
2 * ((int) r1[i + 1] - (int) r1[i - 1]) + \
(int) r2[i + 1] - (int) r2[i - 1]; \
d2 = (int) r2[i - 1] + 2 * (int) r2[i] + (int) r2[i + 1] - \
((int) r0[i - 1] + 2 * (int) r0[i] + r0[i + 1]); \
} while (ZERO)
int
img_edge(IPTR im)
{
int i, j;
register pc_t *r0, *r1, *r2;
register pc_t *g0, *g1, *g2;
register pc_t *b0, *b1, *b2;
long rline;
ci_t **m;
int d1, d2, d;
int gray = IS_GRAY(im);
show_busy("Just a moment...");
if (img_get_rgb(im) < 0)
return -1;
if (!IS_CI(im))
{
/* free the old image raster */
img_free_rastermem(im);
/* new image */
m = get_mat(im->h, im->w, sizeof(ci_t));
}
else
{
m = im->mraster;
}
rline = progress_report("Searching ...", im->h);
memset(m[0], 0, sizeof(ci_t) * im->w);
memset(m[im->h - 1], 0, sizeof(ci_t) * im->w);
for (j = 0; j < im->h; j++)
m[j][0] = m[j][im->w - 1] = 0;
im->cmap->colors = PCMAXV; /* final range */
if (gray)
{
for (j = 1; j < im->h - 1; j++)
{
r0 = get_RM(im)[j - 1];
r1 = get_RM(im)[j];
r2 = get_RM(im)[j + 1];
REPORT(j, rline);
for (i = 1; i < im->w - 1; i++)
{
get_derivative(r0, r1, r2, i);
d = ((float) (Abs(d1) + Abs(d2)) / 2.6);
if (d > im->cmap->colors)
d = im->cmap->colors;
m[j][i] = d;
}
}
}
else
{
for (j = 1; j < im->h - 1; j++)
{
r0 = get_RM(im)[j - 1];
r1 = get_RM(im)[j];
r2 = get_RM(im)[j + 1];
g0 = get_GM(im)[j - 1];
g1 = get_GM(im)[j];
g2 = get_GM(im)[j + 1];
b0 = get_BM(im)[j - 1];
b1 = get_BM(im)[j];
b2 = get_BM(im)[j + 1];
REPORT(j, rline);
for (i = 1; i < im->w - 1; i++)
{
get_derivative(r0, r1, r2, i);
d = Abs(d1) + Abs(d2);
get_derivative(g0, g1, g2, i);
d += Abs(d1) + Abs(d2);
get_derivative(b0, b1, b2, i);
d += Abs(d1) + Abs(d2);
d /= 8;
if (d > im->cmap->colors)
d = im->cmap->colors;
m[j][i] = d;
}
}
}
img_free_rgbmem(im);
/*
* probably should get the map that matches the original image closest in
* terms warmth
*/
get_bit_defmap(!gray, im->cmap->ct[0],
im->cmap->ct[1],
im->cmap->ct[2],
im->cmap->colors);
remove_progress_report();
return fill_image_struct(im, m, im->h, im->w, gray ? T_GMAP : T_CMAP);
}
/********************************************************************
* histogram equalization with thresholding(not written yet)
********************************************************************/
#include "lookup.h"
int
img_enhance(IPTR im)
{
register pc_t *r, *g, *b;
register int row, col, i;
double fact;
long rlines, total = im->h * im->w;
register unsigned int *hist = im->cmap->p_h.hist[3];
static long sum[PCMAX];
if (img_get_rgb(im) < 0)
return -1;
img_get_histogram(im);
memset(sum, 0, sizeof(sum));
sum[0] = hist[0];
for (i = 1; i < PCMAX; i++)
sum[i] = sum[i - 1] + hist[i];
/* get the mapping function */
fact = (double) (PCMAX - 1) / total;
for (i = 0; i < PCMAX; i++)
sum[i] = (fact * sum[i]);
rlines = progress_report("Enhancing Contrast ...", im->h);
r = get_RM(im)[0];
g = get_GM(im)[0];
b = get_BM(im)[0];
for (row = 0; row < im->h; row++)
{
REPORT(row, rlines);
for (col = 0; col < im->w; col++, r++, g++, b++)
{
*r = sum[*r];
*g = sum[*g];
*b = sum[*b];
}
}
remove_progress_report();
row = img_rgb_to_cpack(im);
img_free_rgbmem(im);
return row;
}
/*********************************************************************
* normalize a image. Similar to enhance, the normal curve will based
* on gray levels and comman to all.
**********************************************************************/
int
img_norm(IPTR im)
{
register pc_t *r, *g, *b, *rs;
pc_t lut[PCMAX + 1];
static int crit = 10, hi;
int i;
if (getint("Enter threshold", &crit, 0, (PCMAX - 1) / 2, 0) < 0)
return -1;
if (img_get_rgb(im) < 0)
return -1;
hi = PCMAX - 1 - crit;
for (i = 0; i < crit; i++)
lut[i] = 0;
for (i = hi; i < PCMAX; i++)
lut[i] = PCMAX - 1;
for (i = crit; i < hi; i++)
lut[i] = (hi * (i - crit) / (float) (hi - crit));
r = get_RM(im)[0];
g = get_GM(im)[0];
b = get_BM(im)[0];
for (rs = r + im->w * im->h; r < rs; r++, g++, b++)
{
*r = lut[*r];
*g = lut[*g];
*b = lut[*b];
}
hi = img_rgb_to_cpack(im);
img_free_rgbmem(im);
return hi;
}
/***************************************************************
* Show (RGB seperated) histograms of the current image.
* Note the histogram has two extra entries for max and ave
***************************************************************/
int
img_show_histogram(IPTR im)
{
int i, j;
Hist_t *hist;
float x[PCMAX], y[4][PCMAX], *yy[4], totalpix;
char ave[4][50];
if (!image_ready(im, "Histogram"))
return -1;
if (img_get_histogram(im) < 0)
return -1;
set_chart4_help(HELP_HIST);
hist = im->cmap->p_h.hist;
totalpix = im->h * im->w;
for (i = 0; i < PCMAX; i++)
{
x[i] = i;
for (j = 0; j < 4; j++)
y[j][i] = (double) hist[j][i] / totalpix;
}
for (i = 0; i < 4; i++)
{
yy[i] = y[i];
sprintf(ave[i], "Ave=%u Peak=%.1g at %u",
(unsigned) hist[i][PCMAX + 1], (double) hist[i][PCMAX] / totalpix,
hist[i][PCMAX + 2]);
}
set_chart4_text(ave[0], ave[1], ave[2], ave[3]);
show_chart4("Histogram", x, yy, PCMAX);
return 0;
}
/* ARGSUSED */
int
img_fft(IPTR im)
{
to_be_written("FFT");
return 0;
}
/******************************************************************
* Compute pixel-pixel difference of two identically sized images.
* And show the resulting difference as an image
*******************************************************************/
int
img_diff(IPTR im1)
{
IPTR im2;
register rgba_t *ras2, *ras, *rs;
register int r1, r2, r, g1, g2, g, b1, b2, b;
const char *fname = getfilename("Reference Image", ".", "*.*", "", 1);
if (!fname || !*fname)
return -1;
if ((im2 = open_image(fname)))
{
close_image(im2);
/* currently JPEG desc fakes image size to 1 */
if ((im2->w != 1) && (im2->w != im1->w) && im1->h != im2->h)
{
Bark("IDiff", "%s & %s are of different size", fname,
im1->ifile);
free_image(im2);
return -1;
}
free_image(im2);
im2 = load_image(fname);
}
if (!im2 || !im2->ok)
return -1;
if (im1->w != im2->w || im1->h != im2->h)
{
Bark("IDiff", "%s and %s are of different size", fname, im1->ifile);
free_image(im2);
return -1;
}
strcpy(im1->ifile, "Idiff");
/* both image should be in RGB */
if (img_convert_type(im1, T_RGBA) < 0 ||
img_convert_type(im2, T_RGBA))
{
Bark("IDiff", "type conversion failed");
free_image(im2);
return -1;
}
/* get the difference: put the results into IM1 */
ras = ((rgba_t **) im1->mraster)[0];
ras2 = ((rgba_t **) im2->mraster)[0];
for (rs = ras + im1->w * im1->h; ras < rs; ras++, ras2++)
{
Unpack(*ras, r1, g1, b1);
Unpack(*ras2, r2, g2, b2);
r = Abs(r1 - r2);
g = Abs(g1 - g2);
b = Abs(b1 - b2);
*ras = Pack(r, g, b);
}
free_image(im2);
return 0;
}
/******************************************************************
* Median filtering on a 3x3 grid
***************************************************************{**/
/****************************************************************
* Given RGB matrices and the convolution kernel, do the
* convolution.
****************************************************************/
static int med3x3(pc_t *, pc_t *, pc_t *);
static int
rgb_median(register pc_t **rm, register pc_t **gm, register pc_t **bm,
int nrow, int ncol, int cm_row, int cm_col)
{
register int row, col;
long rlines;
rlines = progress_report("Med3x3", nrow - cm_row / 2);
if (cm_col == 3 && cm_row == 3)
{
nrow--;
ncol--;
for (row = 1; row < nrow; row++)
{
REPORT(row, rlines);
for (col = 1; col < ncol; col++)
{
rm[row][col] = med3x3(rm[row - 1] + col - 1,
rm[row] + col - 1,
rm[row + 1] + col - 1);
gm[row][col] = med3x3(gm[row - 1] + col - 1,
gm[row] + col - 1,
gm[row + 1] + col - 1);
bm[row][col] = med3x3(bm[row - 1] + col - 1,
bm[row] + col - 1,
bm[row + 1] + col - 1);
}
}
}
remove_progress_report();
return 0;
}
static int
gray_median(register pc_t **rm, register pc_t **gm, register pc_t **bm,
int nrow, int ncol, int cm_row, int cm_col)
{
register int row, col;
long rlines;
rlines = progress_report("Med3x3", nrow - cm_row / 2);
if (cm_col == 3 && cm_row == 3)
{
nrow--;
ncol--;
for (row = 1; row < nrow; row++)
{
REPORT(row, rlines);
for (col = 1; col < ncol; col++)
{
rm[row][col] = med3x3(rm[row - 1] + col - 1,
rm[row] + col - 1,
rm[row + 1] + col - 1);
}
}
}
remove_progress_report();
return 0;
}
/**********************************************************************
* Perform median filtering on the part of the image, bounded by mr
**********************************************************************/
int
img_median_filter(IPTR im, Rect_t * mr)
{
rgba_t **subi; /* subimage */
pc_t **rm, **bm, **gm;
int status;
int total = mr->w * mr->h;
/* must be done to a region larger 3x3 */
if (mr->w < 3 || mr->h < 3)
return -1;
/* we can only do this on RGB images */
if (!IS_CPACK(im))
{
img_convert_type(im, IS_GRAY(im) ? T_GRAY : T_RGBA);
/* reset display mode and show image */
im->io->display(im, -1, 0);
}
/* check if on the entire image */
if (!equal_rect(mr, img_rect(im)))
{
if (!(subi = get_subimage(im, mr->x, mr->y, mr->w, mr->h)))
return -1;
rm = get_mat(mr->h, mr->w, sizeof(**rm));
gm = get_mat(mr->h, mr->w, sizeof(**gm));
bm = get_mat(mr->h, mr->w, sizeof(**bm));
/* unpack RGBA to R, G, B. */
unpack_mat(subi[0], rm[0], gm[0], bm[0], total);
}
else
{
img_get_rgb(im);
subi = im->mraster;
rm = get_RM(im);
gm = get_GM(im);
bm = get_BM(im);
}
/* now do the stuff */
status = (IS_GRAY(im) ? gray_median : rgb_median)
(rm, gm, bm, mr->h, mr->w, 3, 3);
/* pack_mat detects gray & rgba by checking if R == G */
pack_mat(subi[0], rm[0], IS_GRAY(im) ? rm[0] : gm[0], bm[0], total);
/* change the image in core */
put_subimage(im, subi, mr, 1);
/* change the image in framebuffer */
set_current_window(win_id);
Rectwrite(mr->x, mr->y, mr->x + mr->w - 1, mr->y + mr->h - 1, subi[0]);
/* clean up */
if (rm != get_RM(im))
{
free_mat(rm);
free_mat(gm);
free_mat(bm);
free_mat(subi);
}
else
{
img_free_rgbmem(im);
}
return status;
}
/*
* Median Finding on a 3-by-3 Grid
* by Alan Paeth
* from "Graphics Gems", Academic Press, 1990
*/
#define s2(a,b) if ((t = b-a) < 0) {a += t; b -= t;}
#define mn3(a,b,c) s2(a,b); s2(a,c)
#define mx3(a,b,c) s2(b,c); s2(a,c)
#define mnmx3(a,b,c) mx3(a,b,c); s2(a,b)
#define mnmx4(a,b,c,d) s2(a,b); s2(c,d); s2(a,c); s2(b,d)
#define mnmx5(a,b,c,d,e) s2(a,b); s2(c,d); mn3(a,c,e); mx3(b,d,e)
#define mnmx6(a,b,c,d,e,f) s2(a,d); s2(b,e); s2(c,f);\
mn3(a,b,c); mx3(d,e,f)
/*
* Find median on a 3x3 input box of integers.
* b1, b2, b3 are pointers to the left-hand edge of three
* parallel scan-lines to form a 3x3 spatial median.
* Rewriting b2 and b3 as b1 yields code which forms median
* on input presented as a linear array of nine elements.
*/
static int
med3x3(pc_t *b1, pc_t *b2, pc_t *b3)
{
register int r1, r2, r3, r4, r5, r6;
register int t;
r1 = *b1++;
r2 = *b1++;
r3 = *b1++;
r4 = *b2++;
r5 = *b2++;
r6 = *b2++;
mnmx6(r1, r2, r3, r4, r5, r6);
r1 = *b3++;
mnmx5(r1, r2, r3, r4, r5);
r1 = *b3++;
mnmx4(r1, r2, r3, r4);
r1 = *b3++;
mnmx3(r1, r2, r3);
return (r2);
}
/***********************************************************************
* End of Median filter
*********************************************************************}*/